In [1]:
using Revise
using ShadyPolytopes
using Polyhedra
using Printf
using LaTeXStrings
using SCS
using Latexify
In [2]:
using WGLMakie, Bonito
Page(offline=true) 
WGLMakie.activate!()

Inscribed and circumscribed sphere of $\mathcal{C}$¶

In [3]:
function euclidean_norm(x)
    return sqrt(sum(t^2 for t in x))
end

function plot_inner_and_outer_sphere(cscb; color=:cornflowerblue)
    vertices = all_vertices(cscb)
    r = minimum(1/euclidean_norm(w) for w in cscb.normals)
    R = maximum(euclidean_norm(v) for v in vertices)
    display(@sprintf "R = %.3f" R)
    display(@sprintf "r = %.3f" r)
    display(@sprintf "R/r = %.3f" R/r)
    display("spectral norm bound squared: $(determine_squared_spectralnorm_bound(cscb))")
    display(@sprintf "spectral norm bound = %.3f" sqrt(determine_squared_spectralnorm_bound(cscb)))
    p = polyhedron(convexhull(convert(Vector{Vector{Float64}},vertices)...))
    n = Polyhedra.Mesh(p)
    fig,ax,P = mesh(n, transparency=true, alpha=0.7, color=color); 
    wireframe!(ax, n, color=:black, transparency=true, alpha=0.3, linewidth=1);
    outer_sphere = mesh!(ax,Sphere(Point3f(0),R), transparency=true, alpha=0.2, color=:grey)
    inner_sphere = mesh!(ax,Sphere(Point3f(0),r), transparency=true, alpha=0.3, color=:black)
    return fig
end
display(L"Analyzing $\mathcal{I}$")
display(ShadyPolytopes.perturbed_icosahedron(-3//5,-1//5,1//10))

fig = plot_inner_and_outer_sphere(optimal_icosahedron);
Analyzing $\mathcal{I}$
CSCP{Rational{BigInt}}(Vector{Rational{BigInt}}[[1, -3//5, 1//10], [1, -1//5, 1//10], [1//10, 1, -3//5], [1//10, 1, -1//5], [-3//5, 1//10, 1], [-1//5, 1//10, 1]], Vector{Rational{BigInt}}[[390//1069, -1250//1069, -710//1069], [710//1069, -390//1069, 1250//1069], [20//53, -55//53, 0], [-710//1069, 390//1069, -1250//1069], [-15//17, 0, -20//17], [-10//9, -10//9, -10//9], [-55//53, 0, 20//53], [-1250//1069, -710//1069, 390//1069], [0, -20//53, 55//53], [-20//17, -15//17, 0], [0, 20//17, 15//17], [-390//1069, 1250//1069, 710//1069], [-20//53, 55//53, 0], [0, -20//17, -15//17], [0, 20//53, -55//53], [1250//1069, 710//1069, -390//1069], [10//9, 10//9, 10//9], [20//17, 15//17, 0], [15//17, 0, 20//17], [55//53, 0, -20//53]])
"R = 1.170"
"r = 0.520"
"R/r = 2.253"
"spectral norm bound squared: 137//27"
"spectral norm bound = 2.253"
In [4]:
fig
Out[4]:

Transformation into approximate John position¶

In [5]:
M = approximate_john_position(optimal_icosahedron)
display(latexify(M))
-------------------------------------------------------------
           Clarabel.jl v0.10.0  -  Clever Acronym              
                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

chordal decomposition:
compact format = on, dual completion = oncronym              
                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

chordal decomposition:

  merge method = clique_graph
  PSD cones initial             = 3
  PSD cones decomposable        = 1
  PSD cones after decomposition = 5
  PSD cones after merges        = 5

problem:
  variables     = 22
  constraints   = 62
  nnz(P)        = 0
  nnz(A)        = 238
  cones (total) = 9
    : Nonnegative = 2,  numel = (1,20)
    : Power       = 2,  numel = (3,3)
    : PSDTriangle = 5,  numel = (6,3,6,10,10)

settings:
  linear algebra: direct / qdldl, precision: Float64
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32ym              
                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

chordal decomposition:

  merge method = clique_graph
  PSD cones initial             = 3
  PSD cones decomposable        = 1
  PSD cones after decomposition = 5
  PSD cones after merges        = 5

problem:
  variables     = 22
  constraints   = 62
  nnz(P)        = 0
  nnz(A)        = 238
  cones (total) = 9
    : Nonnegative = 2,  numel = (1,20)
    : Power       = 2,  numel = (3,3)
    : PSDTriangle = 5,  numel = (6,3,6,10,10)

settings:
  linear algebra: direct / qdldl, precision: Float64
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,

  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,       
                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

chordal decomposition:

  merge method = clique_graph
  PSD cones initial             = 3
  PSD cones decomposable        = 1
  PSD cones after decomposition = 5
  PSD cones after merges        = 5

problem:
  variables     = 22
  constraints   = 62
  nnz(P)        = 0
  nnz(A)        = 238
  cones (total) = 9
    : Nonnegative = 2,  numel = (1,20)
    : Power       = 2,  numel = (3,3)
    : PSDTriangle = 5,  numel = (6,3,6,10,10)

settings:
  linear algebra: direct / qdldl, precision: Float64
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,

  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07

               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0   0.0000e+00  -2.4334e+00  2.43e+00  9.34e-01  4.65e+00  1.00e+00  1.00e+00   ------   l Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

chordal decomposition:

  merge method = clique_graph
  PSD cones initial             = 3
  PSD cones decomposable        = 1
  PSD cones after decomposition = 5
  PSD cones after merges        = 5

problem:
  variables     = 22
  constraints   = 62
  nnz(P)        = 0
  nnz(A)        = 238
  cones (total) = 9
    : Nonnegative = 2,  numel = (1,20)
    : Power       = 2,  numel = (3,3)
    : PSDTriangle = 5,  numel = (6,3,6,10,10)

settings:
  linear algebra: direct / qdldl, precision: Float64
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,

  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07

               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0  
  1  -2.2412e-01  -5.5373e-01  3.30e-01  1.30e-01  9.48e-01  5.21e-02  1.59e-01  8.61e-01  l Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

chordal decomposition:

  merge method = clique_graph
  PSD cones initial             = 3
  PSD cones decomposable        = 1
  PSD cones after decomposition = 5
  PSD cones after merges        = 5

problem:
  variables     = 22
  constraints   = 62
  nnz(P)        = 0
  nnz(A)        = 238
  cones (total) = 9
    : Nonnegative = 2,  numel = (1,20)
    : Power       = 2,  numel = (3,3)
    : PSDTriangle = 5,  numel = (6,3,6,10,10)

settings:
  linear algebra: direct / qdldl, precision: Float64
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,

  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07

               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0  
  1  
  2  -2.6201e-01  -3.7724e-01  1.15e-01  5.11e-02  3.85e-01  2.00e-02  6.35e-02  6.34e-01  
  3  -4.1895e-01  -4.4690e-01  2.79e-02  1.49e-02  8.79e-02  4.82e-03  1.61e-02  7.92e-01  
  4  -4.8033e-01  -4.8614e-01  5.81e-03  3.53e-03  1.76e-02  9.74e-04  3.54e-03  7.92e-01  
  5  -4.8919e-01  -4.8948e-01  2.86e-04  1.76e-04  8.66e-04  5.00e-05  1.77e-04  9.59e-01  
  6  -4.8953e-01  -4.8954e-01  7.90e-06  4.89e-06  2.40e-05  1.42e-06  4.92e-06  9.80e-01  
  7  -4.8954e-01  -4.8954e-01  2.18e-07  1.36e-07  6.66e-07  4.00e-08  1.36e-07  9.80e-01  
  8  -4.8954e-01  -4.8954e-01  6.05e-09  3.76e-09  1.85e-08  1.12e-09  3.78e-09  9.80e-01  
  9  -4.8954e-01  -4.8954e-01  1.30e-09  8.08e-10  3.96e-09  2.41e-10  8.12e-10  7.92e-01  
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time =  8.11s
\begin{equation} \left[ \begin{array}{ccc} \frac{11}{10} & \frac{11}{10} & \frac{11}{10} \\ 0 & \frac{-4}{5} & \frac{4}{5} \\ 1 & \frac{-1}{2} & \frac{-2}{5} \\ \end{array} \right] \end{equation}
In [6]:
display(L"Analyzing $\mathcal{J} = M \mathcal{I}$")
display([M*x for x in optimal_icosahedron.positive_vertices])
display(optimal_icosahedron_john_position.positive_vertices)

fig = plot_inner_and_outer_sphere(optimal_icosahedron_john_position; color=:yellowgreen)
display(fig)
display("1.55 > spectral_norm_bound: ")
display((155//100)^2 > determine_squared_spectralnorm_bound(optimal_icosahedron_john_position))
Analyzing $\mathcal{J} = M \mathcal{I}$
6-element Vector{Vector{Rational{BigInt}}}:
 [11//20, 14//25, 63//50]
 [99//100, 6//25, 53//50]
 [11//20, -32//25, -4//25]
 [99//100, -24//25, -8//25]
 [11//20, 18//25, -21//20]
 [99//100, 18//25, -13//20]
6-element Vector{Vector{Rational{BigInt}}}:
 [11//20, 49//100, 117//100]
 [99//100, 57//100, 81//100]
 [11//20, 81//100, -51//50]
 [99//100, 9//20, -9//10]
 [11//20, -13//10, -3//20]
 [99//100, -51//50, 9//100]
"R = 1.424"
"r = 0.923"
"R/r = 1.543"
"spectral norm bound squared: 370280747942//155558341125"
"spectral norm bound = 1.543"
"1.55 > spectral_norm_bound: "
true

Determining a projection with small norm¶

In [17]:
function makie_polyhedron(vertices)
    vertices = convert(Vector{Vector{Float64}},vertices)
    p = polyhedron(convexhull(vertices...))
    return Polyhedra.Mesh(p)
end

function plot_projection(cscb; color=:cornflowerblue)
    vertices = all_vertices(cscb)
    n = makie_polyhedron(vertices)
    fig,ax,P = mesh(n, transparency=true, alpha=0.7, color=color); 
    wireframe!(ax, n, color=:black, transparency=true, alpha=0.3, linewidth=1);
    proj, shadiness_const = ShadyPolytopes.find_shadiness_constant_numerically(cscb, false)
    approx_proj = ShadyPolytopes.approximate_projection(proj,prec=10^4)
    image = makie_polyhedron([1.2*approx_proj*v for v in vertices])
    mesh!(ax, image, transparency=true, alpha=0.9, color=:darkred);
    kernel = proj*[1,1,0]-[1,1,0]
    lines!(ax, [Point3f(kernel),Point3f(-kernel)],  transparency=true, alpha=0.5,color=:darkred);
    display(@sprintf "s_2 ≈ %.6f" shadiness_const)
    display(@sprintf "s_2 ≤ %.6f" ShadyPolytopes.operator_norm(cscb, approx_proj))
    display("Projection: ")
    display(latexify(approx_proj))
    op_norm = ShadyPolytopes.operator_norm(cscb, approx_proj)
    display("Operator norm of projection = $op_norm ≈ "*(@sprintf "%.5f" op_norm))
    return fig, approx_proj
end

fig, approx_proj = plot_projection(optimal_icosahedron)
fig
presolving:
(round 1, fast)       0 del vars, 1 del conss, 0 add conss, 1 chg bounds, 0 chg sides, 0 chg coeffs, 1 upgd conss, 0 impls, 0 clqs
   (0.0s) symmetry computation started: requiring (bin +, int +, cont +), (fixed: bin -, int -, cont -)
   (0.0s) symmetry computation finished: 1 generators found (max: 1500, log10 of symmetry group size: 0.0) (symcode time: 0.00)
dynamic symmetry handling statistics:
   orbitopal reduction:       no components
   orbital reduction:         no components
   lexicographic reduction:   no permutations
handled 1 out of 1 symmetry components
(round 2, exhaustive) 0 del vars, 1 del conss, 2 add conss, 1 chg bounds, 0 chg sides, 0 chg coeffs, 1 upgd conss, 0 impls, 0 clqs
presolving (3 rounds: 3 fast, 2 medium, 2 exhaustive):
 0 deleted vars, 1 deleted constraints, 0 added constraints, 1 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
 0 implications, 0 cliques
presolved problem has 10 variables (0 bin, 0 int, 0 impl, 10 cont) and 132 constraints
    123 constraints of type <linear>
      9 constraints of type <nonlinear>
Presolving Time: 0.00

 time | node  | left  |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr|  dualbound   | primalbound  |  gap   | compl. 
  0.0s|     1 |     0 |     8 |     - |  1826k |   0 |  43 | 132 | 136 |   0 |  0 |   0 |   0 | 1.000000e+00 |      --      |    Inf | unknown
L 0.2s|     1 |     0 |     8 |     - |  subnlp|   0 |  43 | 132 | 136 |   0 |  0 |   0 |   0 | 1.000000e+00 | 1.281918e+00 |  28.19%| unknown
  0.2s|     1 |     0 |     9 |     - |  1834k |   0 |  43 | 132 | 137 |   1 |  1 |   0 |   0 | 1.000000e+00 | 1.281918e+00 |  28.19%| unknown
  0.2s|     1 |     0 |     9 |     - |  1834k |   0 |  43 | 132 | 137 |   1 |  1 |   0 |   0 | 1.000000e+00 | 1.281918e+00 |  28.19%| unknown
  0.2s|     1 |     0 |  1791 |     - |  1842k |   0 |  43 | 131 | 137 |   1 |  2 |   0 |   0 | 1.000000e+00 | 1.281918e+00 |  28.19%| unknown
  0.2s|     1 |     0 |  1795 |     - |  1858k |   0 |  43 | 131 | 141 |   5 |  3 |   0 |   0 | 1.000000e+00 | 1.281918e+00 |  28.19%| unknown
  0.2s|     1 |     0 |  1797 |     - |  1858k |   0 |  43 | 131 | 144 |   8 |  4 |   0 |   0 | 1.000000e+00 | 1.281918e+00 |  28.19%| unknown
  0.2s|     1 |     0 |  1806 |     - |  1858k |   0 |  43 | 131 | 149 |  13 |  5 |   0 |   0 | 1.000000e+00 | 1.281918e+00 |  28.19%| unknown
  0.2s|     1 |     0 |  1809 |     - |  1858k |   0 |  43 | 131 | 152 |  16 |  6 |   0 |   0 | 1.000000e+00 | 1.281918e+00 |  28.19%| unknown
  0.2s|     1 |     0 |  1811 |     - |  1867k |   0 |  43 | 131 | 154 |  18 |  7 |   0 |   0 | 1.000000e+00 | 1.281918e+00 |  28.19%| unknown
  0.2s|     1 |     0 |  1814 |     - |  1867k |   0 |  43 | 131 | 156 |  20 |  8 |   0 |   0 | 1.000000e+00 | 1.281918e+00 |  28.19%| unknown
L 0.6s|     1 |     0 |  1814 |     - |  subnlp|   0 |  43 | 131 | 156 |  20 |  9 |   0 |   0 | 1.000000e+00 | 1.024936e+00 |   2.49%| unknown
  0.6s|     1 |     0 |  1814 |     - |  1867k |   0 |  43 | 125 | 156 |  20 |  9 |   0 |   0 | 1.000000e+00 | 1.024936e+00 |   2.49%| unknown
  0.6s|     1 |     0 |  1819 |     - |  1867k |   0 |  43 | 125 | 159 |  23 | 10 |   0 |   0 | 1.000000e+00 | 1.024936e+00 |   2.49%| unknown
  0.6s|     1 |     0 |  1822 |     - |  1867k |   0 |  43 | 125 | 157 |  27 | 11 |   0 |   0 | 1.000000e+00 | 1.024936e+00 |   2.49%| unknown
 time | node  | left  |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr|  dualbound   | primalbound  |  gap   | compl. 
  0.6s|     1 |     2 |  1822 |     - |  1867k |   0 |  43 | 125 | 157 |  27 | 12 |   0 |   0 | 1.000000e+00 | 1.024936e+00 |   2.49%| unknown
L 0.7s|    17 |    14 |  2002 | 121.4 |  subnlp|   7 |  43 | 125 | 181 | 133 |  3 |   0 |   0 | 1.000000e+00 | 1.024575e+00 |   2.46%|   0.57%
  1.2s|   100 |    47 |  2956 |  29.3 |  2457k |  15 |  43 | 125 | 160 | 750 |  2 |   0 |   0 | 1.000000e+00 | 1.024575e+00 |   2.46%|   4.72%
  1.4s|   200 |    69 |  4304 |  21.3 |  3101k |  15 |  43 | 125 | 231 |1542 |  3 |   0 |   0 | 1.000000e+00 | 1.024575e+00 |   2.46%|  13.34%
  1.6s|   300 |   111 |  6334 |  21.0 |  4215k |  15 |  43 | 125 | 216 |2856 |  4 |   0 |   0 | 1.000000e+00 | 1.024575e+00 |   2.46%|  13.86%
  1.8s|   400 |   153 |  8135 |  20.2 |  4626k |  17 |  43 | 125 | 254 |4017 |  3 |   0 |   0 | 1.000000e+00 | 1.024575e+00 |   2.46%|  17.58%
  2.0s|   500 |   161 |  9308 |  18.5 |  4802k |  17 |  43 | 125 | 207 |4709 |  4 |   0 |   0 | 1.000000e+00 | 1.024575e+00 |   2.46%|  24.34%
  2.2s|   600 |   213 | 11123 |  18.5 |  5290k |  17 |  43 | 125 | 186 |5884 |  3 |   0 |   0 | 1.000000e+00 | 1.024575e+00 |   2.46%|  30.98%
  2.4s|   700 |   233 | 12179 |  17.3 |  6666k |  17 |  43 | 125 | 204 |6502 |  2 |   0 |   0 | 1.000000e+00 | 1.024575e+00 |   2.46%|  31.99%
  2.6s|   800 |   271 | 13811 |  17.2 |  7108k |  21 |  43 | 125 | 199 |7471 |  2 |   0 |   0 | 1.000000e+00 | 1.024575e+00 |   2.46%|  32.82%
  2.8s|   900 |   287 | 15173 |  16.8 |  7357k |  21 |  43 | 125 | 188 |8349 |  2 |   0 |   0 | 1.000000e+00 | 1.024575e+00 |   2.46%|  38.05%
d 2.9s|   950 |   159 | 16114 |  16.9 |adaptive|  21 |  43 | 125 | 193 |9030 |  3 |   0 |   0 | 1.000000e+00 | 1.012713e+00 |   1.27%|  39.66%
  3.1s|  1000 |   155 | 16651 |  16.6 |  7616k |  21 |  43 | 125 | 193 |9364 |  2 |   0 |   0 | 1.000000e+00 | 1.012713e+00 |   1.27%|  46.47%
  3.4s|  1100 |   155 | 17850 |  16.2 |  7623k |  21 |  43 | 125 | 217 |  10k|  1 |   0 |   0 | 1.000000e+00 | 1.012713e+00 |   1.27%|  52.29%
  3.6s|  1200 |   149 | 19047 |  15.8 |  7639k |  21 |  43 | 125 | 249 |  10k|  2 |   0 |   0 | 1.000000e+00 | 1.012713e+00 |   1.27%|  64.78%
 time | node  | left  |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr|  dualbound   | primalbound  |  gap   | compl. 
  3.8s|  1300 |   153 | 20578 |  15.8 |  7652k |  21 |  43 | 125 | 291 |  12k|  3 |   0 |   0 | 1.000000e+00 | 1.012713e+00 |   1.27%|  68.96%
  3.9s|  1400 |   105 | 21422 |  15.3 |  7662k |  21 |  43 | 125 |   0 |  12k|  0 |   0 |   0 | 1.008217e+00 | 1.012713e+00 |   0.45%|  88.43%
  4.1s|  1500 |    29 | 21974 |  14.6 |  7688k |  21 |  43 | 125 | 266 |  13k|  1 |   0 |   0 | 1.011336e+00 | 1.012713e+00 |   0.14%|  98.38%
* 4.2s|  1537 |     6 | 22137 |  14.4 |    LP  |  24 |  43 | 122 | 350 |  13k|  2 |   0 |   0 | 1.012657e+00 | 1.012713e+00 |   0.01%|  99.99%

SCIP Status        : problem is solved [optimal solution found]
Solving Time (sec) : 4.25
Solving Nodes      : 1547
Primal Bound       : +1.01271284807885e+00 (5 solutions)
Dual Bound         : +1.01271284807885e+00
Gap                : 0.00 %
"s_2 ≈ 1.012713"
"s_2 ≤ 1.012747"
"Projection: "
\begin{equation} \left[ \begin{array}{ccc} \frac{82602121}{79729122} & \frac{54836807}{79729122} & \frac{-722323}{13288187} \\ \frac{-4217717}{79729122} & \frac{-774259}{79729122} & \frac{1060409}{13288187} \\ \frac{695635}{39864561} & \frac{13277555}{39864561} & \frac{12938397}{13288187} \\ \end{array} \right] \end{equation}
"Operator norm of projection = 14386149522//14205071903 ≈ 1.01275"
Out[17]:
In [ ]: